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Introduction 

For over forty years it has been known that, as a result of buoyancy, the ignition of 
near limit hydrogen-air mixtures can lead to the formation of flame-caps-discrete curved 
flames which are small on the scale of the confinement vessel [1]. More recently it has 
been shown that at low gravity, when buoyancy effects are small, flame-balls can be gen- 
erated [2]. These are stationary spherical structures whose existence appears to require a 
near-limit mixture, a small Lewis number and heat losses from radiation. It is our goal 
to combine computational modeling with existing experimental and theoretical studies 
(NASA) of these structures so that an improved understanding of flammability limits and 
near-lifnit phenomena will occur. The question of flammability limits is of fundamental 
importance and has long been examined. It is of great practical importance to predict, 
from first principles, a limit mixture strength that agrees with experimental values for the 
configuration at hand. Flame-balls provide an excellent configuration in which convective 
losses can be eliminated and the resulting stable solutions are produced from a diffusive, 
reactive and radiative balance. 


Problem Formulation 

Although analytical modeling provides convincing evidence that the key physical ingre- 
dients of flame-balls have been identified, quantitative confirmation can only come from 
detailed numerical simulations. Our goal is to predict theoretically the mass fractions 
of the species and the temperature as functions of the independent coordinate r. Upon 
neglecting viscous effects, body forces, and the diffusion of heat due to concentration gra- 
dients, the equations governing the structure of a steady, spherically symmetric, isobaric 
flame-ball are 
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In these equations r denotes the independent spatial coordinate; T, the temperature; Y k> 
the mass fraction of the kfi 0, species; p, the pressure; p, the mass density; W k , the molecular 
weight of the k species; W, the mean molecular weight of the mixture; R, the universal 
gas constant; A, the thermal conductivity of the mixture; c p , the constant pressure heat 

capacity of the mixture; c pk , the constant pressure heat capacity of the k^ species; u k , the 
molar rate of production of the species per unit volume; h k , the specific enthalpy of 
the J species; q R the radiative heat flux and V k , the diffusion velocity of the fc th species. 

Utilizing the kinetic theory of dilute gas mixtures, we can write the following expression 
for the species diffusion velocities 

V k = - f) D kl VX, -6 k V log(T), (4) 

/=1 

where D kh 1 < k,l < K, are the species diffusion coefficients and 6 k , 1 < k < K, are the 
thermal diffusion coefficients. The species diffusion coefficients are symmetric and satisfy 
the important relations 

£y*i> tl = o, i</< k , £y t 0* = o, ( 5 ) 

k = i t = i 

which guarantee the mass conservation constraint 

= ( 6 ) 

*=1 


Using the recent theory of iterative transport algorithms [3], rigorous kinetic theory 
expressions can be derived for all of the transport coefficients. In particular, the thermal 
conductivity, the diffusion coefficients and the thermal diffusion coefficients are obtained by 
solving constrained singular linear systems. These approximate expressions are accurate 
and computationally much more cost-effective than a direct numerical inversion of the 
associated linear systems. Moreover, the approximate species diffusion coefficients and 
thermal diffusion coefficients automatically satisfy the mass conservation constraint (6). 

Since radiation plays a key role in defining the structure and stability of flame-balls, 
it is necessary to construct a careful estimate of the radiation term q R . We assume that 
for hydrogen-air mixtures the only significant radiating species is H^O (For methane-air 
systems additional terms from CO and CO^ will be included). By utilizing an optically 
thin limit in which self absorption of radiation is neglected, the divergence of the net 
radiative flux is given by 

4-£(r 2 4») = 4 < r S>.B J (r) (7) 

r 2 dr , 

where B,(T) is the Planck function evaluated at the band centers of the contributing 
vibration-rotation or pure rotational bands whose integrated intensities are given by q< 
[4,5]. 
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To complete the specification of the problem, boundary conditions are needed at both 
ends of the computational domain. The boundary conditions in the fresh mixture are given 
by 

T(oo) = Ty, n(oo) = y t/ , k = l,2,...,K. (8) 

Symmetry conditions at the origin are given by 

Jrp jy 

* (0) = 0 ’ *r (0) = 0 ' * = 1 ’ 2 >•••.*• (») 

Finally, we point out that the chemical production rates and transport coefficients are 
evaluated with optimized and highly efficient libraries [3,6]. 

Method of Solution 

The equations in (1-4, 7-9) form a nonlinear two-point boundary value problem. The 
solution method for a single solution has been discussed in detail elsewhere and we outline 
only the essential features here [7,8]. Our goal is to obtain a discrete solution of the gov- 
erning equations on the finite difference mesh M.. By replacing the continuous differential 
operators by finite difference expressions, the problem of finding an analytic solution of the 
governing equations is converted into one of finding an approximation to this solution at 
each point of the mesh M.. We seek the solution U m of the nonlinear system of difference 
equations 

F(U) = 0. (10) 

Assuming we can obtain an initial solution estimate U° that is sufficiently “close” to U*. 
the system of equations in (10) can be solved by Newton’s method. We write 

J(U k ) (U k+1 - U k ) = -\ k F(U k ), k = 0,1,..., (11) 

where U k denotes the k^ 1 solution iterate, A t the k^ 1 damping parameter (0 < A < 1) and 
J{U k ) = dF(U k )/dU the Jacobian matrix. A system of linear block tridiagonal equations 
must be solved at each iteration for corrections to the previous solution vector. As we 
found in the solution of burner-stabilized and freely propagating premixed laminar flames, 
the cost of forming (we use a numerical Jacobian) and factoring the Jacobian matrix can 
be a significant part of the cost of the total calculation. In such problems we apply a 
modified Newton method in which the Jacobian is re-evaluated periodically. 

The solution of combustion problems, such as the flame-ball system, requires that the 
computational mesh be determined adaptively. We equidistribute the mesh M. on the in- 
terval [0, oo] with respect to the non-negative function W and the constant C. Specifically, 
we form the mesh by employing a weight function that equidistributes the difference in 
the components of the discrete solution and its gradient between adjacent mesh points. 

In addition to the generation of a single solution for these systems, we have studied 
the dependence of these flames (the flammability limits) on the equivalence ratio <j>. While 
the solution method discussed above can be utilized to generate a single solution, it is 
not effective for the computation of solutions as one approaches turning points of the 
system. To be able to study flames in the neighborhood of extinction, we must apply a 
phase-space, pseudo-arclength continuation method with Newton-like iterations and global 
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adaptive gridding [9-11]. The set of governing equations for the flame-ball configuration 
reduces to a system of the form 

r{XA) = 0 , ( 12 ) 

where X = (T, Y x , . . . , Yk) is the solution vector. The solutions (X, <f>) in (12) form a 
one-dimensional manifold which, as a result of the presence of turning points, cannot be 
parameterized in the form ( X(<f > ), 4>). The upper part of the manifold denotes the stable 
solutions and the lower part the unstable ones assuming there are no Hopf bifurcations. 

To generate this solution set, ( X,<j> ) is reparameterized into (X(s),<f>($)) where s is a 
new independent parameter and <f> becomes an eigenvalue. The system in (12) can now be 
written 

JF(*(s),^(s))=0, (13) 

and the dependence of s on the augmented solution vector ( X , <f>) is specified by an extra 
scalar equation 

M(X(s),t(s),s)=0, (14) 

which is chosen such that s approximates the arclength of the solution branch in a given 
phase space. Rather than solving the coupled system in (13) and (14), we replace <f> by 
a f un ction <f> of z and we let the unknown 2. = ( X T , <f>) T be the solution of a three-point 
limit value problem 

’*{*,*) 

H(Z,s)= g =0, (15) 

s). 

where z is a given point in [0, oo]. The system in (15) is solved by combining a first- 
order Euler predictor and a corrector step involving Newton-like iterations and adaptive 
gridding. 

Numerical Results 

During the first eight months of the project our goals have included verifying the 
sensitivity of the flame ball radius and the value of the lean and rich flammability limits to 
1) the hydrogen-air chemistry, 2) the transport model and 3) the radiation approximation. 
Otar previous hydrogen-air flame ball computations [12] employed a Curtiss-Hirschfelder 
transport approximation with a radiation model in which interpolated emissivity data was 
used to construct an approximation for the divergence of the radiative heat flux. 

Utilizing a modified hydrogen-air reaction mechanism [13] with the transport model 
employed in [12], we computed flame balls as a function of the equivalence ratio <j>. Tem- 
perature and radii results are compared with our original computations in Figures la and 
lc. While the temperature and the size of the flame ball at the rich limit show only minor 
changes, there is a modest increase of the rich flammability limit with the new mechanism. 
Only very small changes occur on the lean side but they move in the direction of increased 
agreement with the experiments [12]. In Figures lb and Id we compare the temperature 
and radii as a function of the equivalence ratio for the new mechanism/ old transport cal- 
culations with calculations made with the new mechanism and the transport and radiation 
models discussed above. Essentially no major differences were found. Nevertheless, due to 
its more accurate formulation, the detailed transport and band radiation models will be 
utilized in all future flame ball calculations. 
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Figure 1 (a) A comparison of the temperature versus the equivalence ratio for hydrogen- 
air flame balls utilizing the old mechanism and old transport (dash) with computations 
utilizing the new mechanism and old transport (solid); (b) a comparison of the temper- 
ature versus the equivalence ratio for the new mechanism and old transport (dash) with 
computations utilizing the new mechanism, transport and radiation models (solid); (c) a 
comparison of the flame ball radius versus the equivalence ratio for the old mechanism 
and old transport (dash) with computations utilizing the new mechanism and old trans- 
port (solid); (d) a comparison of the flame ball radius versus the equivalence ratio for the 
new mechanism and old transport (dash) with computations utilizing the new mechanism 
transport and radiation models (solid). 
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